# -*- coding: utf-8 -*-
"""
Created on Mon Sep 17 10:27:11 2012

@author: xudi
"""

import numpy as np
import pylab as pl

size=6000
maxtime=100000
N=10
C=3e8
lambc=7.8e-7
dz=lambc/1000
dt=dz/(2*C)
e0=8.85e-12
sampling_rate=1/dt

tw=5e-16# only one can change 5e-16
t0=2e-15
t = np.arange(0.,maxtime*dt,dt)
th=t-0.5*dt

wc=2.417e15
w0=100*wc
u0=1.26e-6
Tao=wc
value1=(1-(Tao*dt/2))/(1+(Tao*dt/2))
value2=dt/(1+(Tao*dt/2))

Hy=np.zeros(size)
Ex=np.zeros(size)
R=np.zeros(maxtime)
T=np.zeros(maxtime)
Jx=np.zeros(size)
Px=np.zeros(size)    

def medium(Num,SIZE):
    xl1=1
    xl2=3
    e1=np.zeros(177)
    e2=np.zeros(125)
    e1[:177]=xl1*e0+e0
    e2[:125]=xl2*e0+e0    
    M=np.zeros(SIZE)
    e=np.append(e1,e2)
    mm=0
    M[:1000]=e0
    mm+=1000
    for mm in range(1000,SIZE):
        if (mm-1000)%302==0:
            M[mm:mm+302]=e
            mm+=302
            if mm>1000+302*Num/2-1:
                M[mm:mm+177]=xl1*e0+e0  #carefull
                mm+=177
                M[mm:SIZE]=e0
                break
        mm-=1
    #print mm
    return M

Mdm=medium(N,size)   

def source(T):
    ft=np.exp(-(T-t0)**2/tw**2)
    return ft


Xs=source(t)
Ys=source(th)
    
for qtime in range(0,maxtime):
    
    UHy0=Hy[0]
    UHy1=Hy[1]
    UHy2=Hy[size-1]
    UHy3=Hy[size-2]
    
    Hy[1:size]=Hy[1:size]-(dt/(u0*dz))*(Ex[1:size]-Ex[:size-1])
    
    Hy[0]=UHy1+((C*dt-dz)/(C*dt+dz))*(Hy[1]-UHy0)
    Hy[size-1]=UHy3+((C*dt-dz)/(C*dt+dz))*(Hy[size-2]-UHy2)
    
    Hy[700]+=(dt/(u0*dz))*Xs[qtime]
    
    Jx[:size]=value1*Jx[:size]+value2*(Mdm[:size]*(wc**2)*Ex[:size]-(w0**2)*Px[:size])
  
    #UEx0=Ex[0]
    #UEx1=Ex[1]
    Ex[:size-1]=Ex[:size-1]-(dt/Mdm[:size-1])*((Hy[1:size]-Hy[:size-1])/dz+Jx[:size-1])
    
    #Ex[0]=UEx1+((C*dt-dz)/(C*dt+dz))*(Ex[1]-UEx0)
    Ex[700]+=((dt/(e0*dz))*(e0/u0)**0.5)*Ys[qtime]
    #Ex[1000]+=Xs[qtime]
    Px[:size]=Px[:size]+dt*Jx[:size]
    
    R[qtime]=Ex[100]
    T[qtime]=Ex[5900]


fft1=np.fft.fft(T,maxtime)
fft2=np.abs(np.fft.fftshift(fft1))
fft3=np.fft.fft(R,maxtime)
fft4=np.abs(np.fft.fftshift(fft3))
fft5=np.fft.fft(Xs,maxtime)
fft6=np.abs(np.fft.fftshift(fft5))
freqs=np.linspace(0,sampling_rate,maxtime)-sampling_rate/2
y1=np.arange(0,maxtime*dz,dz)

pl.figure(figsize=(16,6))
pl.subplot(221) 
pl.plot(freqs,fft6,'r')
pl.xlim(0,2*wc)
pl.xlabel(" frequence (Hz)")
pl.ylabel("Gauss FFT")


pl.subplot(222)
pl.plot(y1,T,'r')
pl.plot(y1,R,'b')
pl.xlabel("Time (s)")
pl.ylabel("Gauss inedx")
pl.legend((u"T ","R "))

pl.subplot(223)
pl.xlim(0,2*wc)
pl.plot(freqs,fft2,'r')
pl.plot(freqs,fft4,'b')
pl.xlabel(" frequence (Hz)")
pl.ylabel("Reflection and Transmission inedx")
pl.legend((u"T (fft)","R (fft)"))

pl.subplot(224)
pl.xlim(0,1e15)
pl.ylim(0,1)
pl.plot(freqs,fft2/fft6,'r')
pl.plot(freqs,fft4/fft6,'b')
pl.xlabel(" frequence (Hz)")
#pl.ylabel("Refractive inedx")
pl.legend((u"T/G","R/G "))

